Computational studies of the glass-forming ability of model bulk metallic glasses 
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Bulk metallic glasses (BMGs) are produced by rapidly thermally quenching supercooled liquid 
metal alloys below the glass transition temperature at rates much faster than the critical cool- 
ing rate Rc below which crystallization occurs. The glass-forming ability of BMGs increases with 
decreasing Rc, and thus good glass-formers possess small values of Rc- We perform molecular dy- 
namics simulations of binary Lennard- Jones (LJ) mixtures to quantify how key parameters, such 
as the stoichiometry, particle size difference, attraction strength, and heat of mixing, influence the 
glass-formability of model BMGs. For binary LJ mixtures, we find that the best glass-forming 
mixtures possess atomic size ratios (small to large) less than 0.92 and stoichiometries near 50:50 
by number. In addition, weaker attractive interactions between the smaller atoms facilitate glass 
formation, whereas negative heats of mixing (in the experimentally relevant regime) do not change 
Rc significantly. These studies represent a first step in the development of computational methods 
for quantitatively predicting glass-formability. 
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INTRODUCTION 



When supercooled liquids are rapidly quenched at 
rates R exceeding a critical value Rc, crystallization is 
avoided, and systems form disordered solids such as bulk 
metallic glasses (BMGs). BMGs possess high mechani- 
cal strength and can be processed so that they display 
plastic [1|, not brittle, response to applied deformations, 
which makes them desirable materials for a variety of in- 
dustrial and engineering applications [2j]. Avoiding crys- 
tallization in pure metals requires enormously large cool- 
ing rates in excess of 10^^ K/s. However, bulk metallic 
glass-forming alloys have been developed for which the 
critical cooling rate is more than nine orders of magni- 
tude lower, in the range 1 < i?c < 10^ K/s. Understand- 
ing the important physical quantities that determine the 
glass-forming ability of multi-component alloys will allow 
us to develop even stronger and less costly bulk metallic 
glasses. 

Prior research suggests that multi-component metallic 
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alloys with Tg/T^ > 0.67 form BMGs, 
are the glass transition and melting temperature, respec- 
tively [3|. In addition, Inoue 2] has emphasized three 
guidelines for enabling BMG formation, rather than crys- 
tallization: 1) atomic size ratios (small relative to large) 
of a < 0.89 for at least two constituents of the alloy; 2) 
large negative heats of mixing [^]; and 3) several atomic 
components. In Fig. [1] we show the distributions of the 
atomic size ratios and heats of mixing for common bi- 
nary and ternary bulk metallic glass- forming alloys [5|. 
For binary systems, the most probable atomic size ratio 
is a ~ 0.8 and heat of mixing is negative and roughly 



6-7% of the average cohesive energy. 

However, beyond these heuristic guidelines, there is no 
quantitative and predictive understanding of the glass- 
forming ability in multi-component alloys. (Note that 
there have been previous measurements of the critical 
cooling rate in binary hard-sphere systems [1, For 
model BMG-forming systems with attractive interac- 
tions, we do not know the dependence of the critical 
cooling rate on the stoichiometry, size ratios, and heats 
of mixing of the constituent atomic species. For example, 
can multi-component systems with large negative heats 
of mixing, but smaller atomic size mismatches possess the 
same glass-forming ability as systems with small negative 
heats of mixing but larger atomic size mismatches? 

We perform molecular dynamics simulations of model 
glass-forming systems, binary Lennard- Jones mixtures 
of spherical particles, to measure the critical cooling 
rate as a function of the size ratio, number fraction, 
and interaction energy of the two particle species. We 
find several important results. First, the critical cooling 
rate decreases exponentially with the particle size ratio, 
Rc ~ exp[— C(l — a)^], where C depends on the number 
fraction of small and large particles. At a given size ra- 
tio a < 1, the minimum critical cooling rate occurs at 
the number fraction corresponding to equal volumes of 
the large and small particles. In addition, we find that 
at fixed number fraction and size ratio, the critical cool- 
ing rate decreases strongly with decreasing cohesive en- 
ergy ratio of the small particles relative to the large ones, 
cbb/^aa- In contrast, variations of the heat of mixing 
of the two species in the experimentally accessible range 
do not affect Rc significantly. Thus, we have quantified 
several design principles for improving glass formation in 
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FIG. 1: (color online) (left) Probability distribution P{(Jbb /o'aa) of atomic size ratios (with asB < ""aa) in binary bulk 
metallic glasses 0|. (middle) Probability distributions of atomic size ratios P{(Jbb / ctaa) (shaded) and P{acc / ctaa) (white) 
(with acc < o"ss < ctaa) in ternary BMGs [5j]. (right) Probability distribution of the heats of mixing AHmix relative to the 
average cohesive energy (eaa + ess )/2 in binary BMGs [l,!!]. 



binary mixtures. 



2. SIMULATION METHODS 

We perform constant number, volume, and temper- 
ature (NVT) molecular dynamics (MD) simulations of 
binary Lennard- Jones (LJ) mixtures of iV = Na + Nb 
spherical particles with the same mass to, but different 
diameters <jaa and asB, in periodic cubic cells with vol- 
ume V = . The particles interact pairwise via the LJ 
potential 
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where i,j e {A,B}, B indicates the smaller particle, 
<Jij — {<yii + ^jj)/'^' unless otherwise specified, and caa 
and tBB represent the cohesive energies for the A and B 
particles, respectively. We quantify the heat of mixing 
using AiJinix = {€AA + ^bb)/'^ - GAB- We employ the 
shifted-force version of the LJ potential (Eq. [IJ so that 
the pair potential and force vanish for separations beyond 
the cutoff distance rcut = 3.5<7ij [8]. Energies, lengths, 
timescales, and temperatures are given in units of saa, 
caa, cTAA^/m/eAA, and eAA/kB, respectively, where the 
Boltzmann constant ks is set to unity. 

We study the glass-forming ability of binary LJ mix- 
tures at fixed packing fraction ~ A^cr^^(l -I- fB[oi^ — 
l))7r/6y = 0.5236 as a function of the number fraction 
Jb = Nb/N, particle size ratio a — (Jbb/<^aa, relative 
cohesive energy tBB/^AA, and heat of mixing Ai/mix. 
We only show results for 0.92 < a < 1 for which solid 
solutions with FCC crystal structures are the equilibrium 
phase . We initialize the systems at high temperature 
To = 2.0, using the Nose- Hoover thermostat [id , [ll|, 
and then thermally quench the systems exponentially, 
T{t) = Toe-^\ from To to T/ = 10"^ at various rates 
R over four orders of magnitude. (In Appendix \^ we 
show that our results are not sensitive to the choice of 
the thermostat and the form of the cooling schedule.) 

Following the thermal quenches to Ty, we characterize 



the structural properties of the system by measuring sev- 
eral quantities: 1) the local and global bond orientational 
order parameters [l2| - [l3 | 
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where 9ij and 



ij are the axial and polar angles between 
each particle i and its neighbors j, YjJ" are spherical har- 
monics of degree 6 and order to, and is the number of 
nearest neighbors of particle i within a cutoff distance of 
1.5(7^ ; 2) local bond orientational order position correla- 
tion function 
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where g{r) — '^i'^j^i^^r — rij) is the radial distribu- 
tion function and g6m(n) = Yl%iYi^{Oij , and 
3) the crystal domain size. These structural quantities 
are averaged over at least 96 independent quenching tra- 
jectories. (In Appendix iBl we compare the results using 
these structural quantities.) We consider system sizes 
from N = 500 to 8788 particles. 



3. RESULTS 

In this section, we characterize the structural proper- 
ties of LJ systems thermally quenched to temperature Tj 
as a function of the cooling rate R. In the right inset of 
the left panel of Fig. [21 we show the distribution P{Qq) 
of the local bond orientational order parameter Qg for 
monodisperse LJ systems with N = 1372 particles. For 
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FIG. 2: (color online) (left) Median local bond orientational order parameter Qg for monodisperse Lennard- Jones (LJ) systems 
following thermal quenches to Tf = 0.01 over a range of cooling rates R for system sizes A'^ = 500, 864, 1372, 2048, 4000, 

and 8788. The critical cooling rate Rc (defined, as discussed in the main text, as the rate at which Qg = 0.43 (dashed line)) 
approaches its large- A'^ limit, 7?^, as a power-law Rc — RT ~ l/N"^ (left inset), (right inset) The probability distribution ^(Qe) 
for monodisperse LJ systems with N = 1372 following quenches to Tf — 0.01 for cooling rates R = 0.02 (o), 0.01 (A), and 
0.005 (□). (right) Fraction of particles that occur in HOP (o), FCC (□), and BCC (A) crystal clusters as a function of 1/A^ 
for monodisperse LJ systems following a quench to Tf at cooling rate R = 10"'' < Rc- At this rate, the median local bond 
orientational order parameter Qg (x) agrees with the value (V) obtained by averaging Qg = 0.575 for FCC and = 0.485 for 
HCP (dashed lines) weighted by the fraction of particles in FCC and HCP clusters in each sample. 



fast cooling rates, e.g. R = 0.02, most of the quenched 
systems are structurally disordered, and P{Q^q) possesses 
a strong peak at small Qg ^ 0.41. In contrast, for slow 
cooling rates, e.g. R — 0.005, most of the quenched sys- 
tems are ordered, and P{Q\) possesses a strong peak at a 
larger value of Qg ~ 0.51. For intermediate cooling rates, 
the distribution P{Q\) becomes strongly bimodal, which 
indicates that the systems possess disordered as well as 
ordered regions. In the main panel of Fig. [2] (left), we 

show the median Qg versus the logarithm of the cooling 

rate R for several system sizes. For each system size, Qg 
first increases modestly with decreasing cooling rate, fol- 
lowed by a rapid increase at intermediate rates, and then 
it plateaus with further decreases. We define the critical 
cooling rate, as the rate at which the median local 
bond orientational order parameter crosses the thresh- 
old value Qq = Qo = 0.43. We chose the threshold Qo 

for several reasons: 1) Qo captures the steep rise in Qg 
with decreasing cooling rate, 2) Qq is in the region of Qg 
between the two peaks in P{Q\) that occur at intermedi- 
ate cooling rates (right inset of left panel of Fig. [5]), and 
3) Qo is a value for which Q\,{R) becomes system size 
independent for intermediate and fast cooling rates. 

Note that the distribution of the global bond orienta- 
tional order parameter P{Qq) also becomes bimodal and 
the median Qg increases rapidly with decreasing cooling 
rate. (See Appendix B.) However, the global bond ori- 
entational order parameter quantifies crystallization of 
the entire system, which is influenced more by the slow 
dynamics of crystal growth, rather than the initial nucle- 




FIG. 3: (color online) Crystalline clusters obtained in 
monodisperse LJ systems with A'^ = 4000 following ther- 
mal quenches to Tf — 0.01 at cooling rates R = 10~^ (left) 
and 10~^ (right). Particles are colored according to whether 
they belong to FCC (cyan), HCP (blue), BCC (red), or non- 
crystalline (white) domains. 



ation of crystalline domains. 

The value of the bond orientational order parameter 
depends on the crystal structure that forms during the 
thermal quenching process. Thus, we employed a crys- 
tal analysis algorithm to identify the crystalline clusters 
(FCC, HCP [151], or BCC) for cooling rates R < Rc- 
For example, Qg ~ 0.575 for an ideal face-centered cu- 
bic (FCC) structure, whereas it is « 0.485 for an ideal 
hexagonal close packed (HCP) structure. This difference 

explains the increase in Qg for R <SC_ Rc as N increases 
in the main panel of Fig. [5] (left). In Fig. [5] (right), we 
show that small systems N < 500 mainly crystallize to 
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FIG. 4: (color online) (left) Critical cooling rate Rc for binary LJ mixtures with A*' = 1372 as a function of the number fraction 
fs for several particle size ratios a = 1.0, 0.97, 0.95, 0.94, 0.93, and 0.92. The solid lines are sixth-order least-square fits to 
the data for Rc- The dashed line connects the number fractions fg = 1/(1 -|- a^) at which the A and B particles occupy the 
same volume, (right) Rc versus (1 — a)^ for binary LJ mixtures with fs = 0.2, 0.4, 0.5, 0.6, and 0.8. The error bars for Rc 
are determined by the cooling rate increment AR = 10~^. The inset shows the coefficient C{fB) of the exponential decay of 
Rc ~ exp[-C(l - af\. 



HCP structures [T^ , while larger systems crystallize pre- 
dominantly to FCC structures. For low cooling rates, the 
median local bond orientational order parameter can 
be obtained by averaging the Qg values for FCC and HCP 
structures weighted by the fraction of particles in FCC 
and HCP clusters in each sample. (See Fig. [2] (right).) 
We show snapshots of the thermally quenched structures 
for monodisperse LJ systems using two cooling rates in 
Fig. [3] with FCC, BCC, HCP, and non-crystalhne regions 
shaded different colors. 

We show the system-size dependence of the critical 
cooling rate Rc for monodisperse LJ systems in the left 
inset to the left panel of Fig. We find that Rc decreases 
with increasing system size and approaches its large- 
hmit, Rf w 0.01, as a power law Rc - Rf ^ l/N'^. It 
is interesting that the approach to i?^ scales as l/N'^, 
which is faster than the 1/N scaling typical for first-order 
transitions. In contrast to hard-sphere systems [13], crys- 
tallization in monodisperse L J systems is more difficult at 
large N. In small monodisperse LJ systems {N < 500), 
the critical nucleus is sufficiently large that it interacts 
with its periodic images , which reduces the inter- 

facial energy of crystal nuclei and enhances the formation 
of single crystals. 

We now focus on binary LJ systems at fixed N = 1372 
and cohesive energy ratio eBs/^AA = 1 and study the 
glass-forming ability as a function of the size ratio a 
and number fraction fs- For a < 1, the smallest 
Rc{oi, Jb) (j-e. best glass-former) is obtained in systems 
with approximately equal numbers of A and B particles, 
« 0.5, as shown in Fig. [4] (left). As a decreases, 
the minimum in Rc{a, fs) deviates from fg w 0.5 and 
follows = 1/(1 -I- a^) for which the A and B par- 
ticles occupy the same volume (reaching fg w 0.56 at 



a — 0.92). As shown in Fig. 2] (right), at each fs, 
Rc decreases exponentially with decreasing size ratio, 
RciaJs) = i?c(l,/s)exp[-C(/s)(l - a)^]. This resuh 
implies that Rc drops from 10^^ to 10~^^-10~^^ for bi- 
nary systems of composition = 0.2-0.8 with size ra- 
tio a = 0.8 (the most common size ratio in binary bulk 
metallic glass formers), which is 9-23 orders slower than 
the Rc aX a — 1. We also note that for a given cooling 
rate R, the glass- forming regime, i.e. the range of num- 
ber fractions for which R > Rc, expands with decreasing 
a. 

For the results presented so far, we set the cohesive 
energy ratio cbb/saa = 1- However, as shown in the 
inset to Fig. [SI the cohesive energy between like species 
is different for the two components for most binary bulk 
metallic glass formers. In Fig.[5l we show Rc as a function 
of ^bb/^aa for binary LJ mixtures with N = 1372 at 
fixed /b = 0.5, (cAA + ^bb)/'2' ~ 1, and heat of mixing 
ATJmix = 0, assuming AiJmix = (eAA+eBs)/2-eAS 0- 
22] and the mixing rules tAB = (eAA+eBB)/2 and (Tab = 
{cTAA + <^bb)/2. We find that the glass-forming ability 
increases {i.e. Rc decreases) as cbb/^aa decreases below 
1 . This result is consistent with the fact that most binary 
glass formers with 0.8 < a < 1 possess ^bbI^aa < 1 [a]- 
(See the inset to Fig.O) 

Inoue's guidelines |2| suggest that a negative heat of 
mixing ATJ^ix < enhances the glass-forming ability of 
BMGs. The rationale is that a negative heat of mix- 
ing makes the mixed and geometrically frustrated state 
energetically favorable compared to the phase separated 
state. Fig. [T] (right) shows that AHj^ix is approximately 
5-10% of the average cohesive energy of the two compo- 
nents, (eaa + ^bb)/'^, for most binary BMGs 0, H, HI]- 
However, we show in Fig. [5] that binary LJ mixtures with 
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FIG. 5: (color online) Critical cooling rate Rc as a function 
of the cohesive energy ratio eBs/^AA for binary LJ mixtures 
with A'^ = 1372, number fraction fs ~ 0.5, and size ratios 
a = 1.0, 0.97, and 0.95. The solid lines indicate results for 
the mixing rule ctas = a = {(jaa + (Jbb)/'^- and €ab = e = 
{(-AA + tBB^I'i- The dashed lines indicate results for positive 
{tAB = 0.9e) (V) and negative heats of mixing (e^s = l.le) 
(A) with GAB = tr and for bond shortening (tab = 0.99a- with 
A/Zmix = (❖). (inset) Cohesive energy ratio tBB/f-AA versus 
the atomic size ratio obb/(taa for common binary metallic 
glass formers [^. 

heats of mixing in the range 2Ai?i„ix/ {f^AA+f-Bs) = ±0.1 
possess the same critical cooling rate Rc as those with 
Ai?niix = over the full range of size ratios studied. 

Why then do most BMGs possess AiJmix < 0? One 
possibility is that negative heats of mixing are corre- 
lated with strong bonding between atomic species, which 
can be modeled as bond shortening {jtab < {o'aa + 
aBB)/'2) l2ji26i]- In Fig. [H we show that only a 1% 
bond shortening, ctas — 0.99((Taa + (^bb)/^, can give 
rise to a finite decrease in the critical cooling rate Rc- 

4. CONCLUSION 

The glass formability of bulk metallic glass-forming al- 
loys can be characterized by the critical cooling rate Rc, 
below which the system possesses crystalline domains. 
The best bulk metallic glasses are those with the lowest 
values for Rc- However, the key parameters that deter- 
mine Rc are not currently known, and thus BMGs are 
mainly developed through a trial and error process. As 
a first step in computational design of BMGs, we per- 
formed molecular dynamics simulations of coarse-grained 
models for BMGs, binary Lennard- Jones mixtures, and 
measured Rc as a function of the number fraction, size ra- 
tio, relative cohesive energy, and heat of mixing of the two 
atomic species. We measured the local bond orientational 
order parameter to quantify the degree of crystallization 
that had occurred in systems during thermal quenches 
from high to low temperature over more than four or- 



ders of magnitude in the cooling rate. It is known that 
weakly polydisperse LJ systems are poor glass-formers; 
we quantified this statement by showing that the criti- 
cal cooling rate decreases exponentially with increasing 
particle size ratio a, Rc ^ exp[— C(l — a)^]. Further, 
at a given size ratio a < 1, the minimum critical cool- 
ing rate occurs at the number fraction corresponding to 
equal volumes of the large and small particles of equal 
mass. In addition, we find that at fixed number fraction 
and size ratio, the critical cooling rate decreases strongly 
with decreasing cohesive energy ratio of the small parti- 
cles relative to the large ones, cbb/^aa- This result may 
explain why most experimentally obtained binary BMGs 
possess ^bb/^aa < 1- In contrast, variations of the heat 
of mixing of the two species in the experimentally ac- 
cessible range (several per cent of the average cohesive 
energy) do not affect Rc for binary LJ mixtures signifi- 
cantly. However, bond shortening of only several percent 
relative to (Tab = {<^aa + (^bb)/'^ [24-26] does give rise 
to significant changes in Rc. Recent experiments have 
suggested that negative heats of mixing are correlated 
with bond-shortening, which may explain why most ex- 
perimentally obtained BMGs possess negative heats of 
mixing. In future studies, we will characterize the glass- 
forming ability and crystallization processes in ternary 
and quaternary LJ mixtures using MD simulations, en- 
ergy minimization, and genetic algorithms. 
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Appendix A: Thermostat and Quenching Protocol 

In this appendix, we provide additional details of the 
molecular dynamics (MD) simulations used to thermally 
quench Lennard-Jones (LJ) systems from high tempera- 
ture liquids to low temperature glasses. The LJ liquids 
were first equilibrated at high temperature Tg = 2.0 using 
constant number N , volume V , and temperature T MD 
simulations, and cooled exponentially T{t) = Toe~^'^ to 
low temperature Tf = 10~^. The temperature was con- 
trolled using the Nose-Hoover thermostat [13, [HI with 
thermal inertia parameter Q = 1, and the equations of 
motion were integrated using a Newton's method tech- 
nique [13 with time step At = 10^^. In Fig. [HI (left), we 
show for monodisperse LJ systems with N = 4000 that 
the dependence of the median local bond orientational 

parameter Qq on rate R is the same for Q = 1 and 10. 

We also investigated the extent to which the thermo- 
stat affects the critical cooling rate, below which the sys- 
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FIG. 6: (left) Median local bond-orientational order parameter q\ versus the cooling rate R for monodisperse LJ systems 
with A'^ — 4000 using the Nose-Hoover thermostat with thermal inertia parameter Q = 1 (□) and 10 (0) in units of ma\j^. 

(middle) Median local bond-orientational order parameter Qg versus R for monodisperse LJ systems with TV — 500 using 
several thermostats: Nose-Hoover (□), Gaussian constraint (0), and ad hoc velocity rescaling A. (right) Median local bond- 
orientational order parameter Qg versus the cooling rate R for monodisperse systems with A'^ = 1372 for a linear thermal 
quenching protocol, T{t) = To — Rt. 
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FIG. 7: (color online) Median global bond orientational order 
parameter Qg for monodisperse LJ systems following thermal 
quenches to Tf = 0.01 over a range of cooling rates R for 
system sizes N = 500, 864, 1372, 2048, 4000, and 8788. (in- 
set) The probability distribution P(Q^) for monodisperse LJ 
systems with A'' = 1372 following quenches to Tf = 0.01 for 
cooling rates R = 0.02 (o), 0.01 (A), and 0.005 (□). 



terns crystallize. In Fig. [HI (center), we show that Qg 
versus R is the same for monodisperse LJ systems with 
N — 500 when the temperature is controlled using the 
Nose-Hoover, Gaussian constraint, and ad hoc velocity 
rescahng thermostats H^. Thus, the choice of the 
thermostat does not influence the measurement of Re- 
We also varied the form of the thermal quenching pro- 
tocol. In Fig. El (right), we show that a linear cooling 
schedule, T{t) = Tq — Rt, gives qualitatively the same 

results for Qq versus R as an exponential temperature 
ramp. 



Appendix B: Characterization of Crystalline Order 

In this Appendix, we describe several metrics (in ad- 
dition to the local bond orientational order parameter 
Qq in Eq. ^ to characterize the degree of crystalline or- 
der of thermally quenched LJ systems. In contrast to 
Qg, the global bond orientational order parameter in 
Eq. [3] quantifies the degree of crystallization over the en- 
tire system. The median global bond orientational order 
parameter Qg versus cooling rate R for monodisperse LJ 
systems for several system sizes is shown in Fig. [T] Qg 
shows a rapid increase near the critical cooling rate Rc as 
found for Qg. However, Rc (defined by a threshold such 
as Qg = 0.3) appears to decrease to zero in the large 
system limit. This trend occurs because it takes an in- 
creasing amount of time (and thus slower cooling rates) 
for crystal nuclei to grow and for the system to reach the 
same Qq as that obtained in smaller systems. 

In Fig. |8] (left), we show the local bond orientational 
order correlation function (Eq. \^ for monodisperse LJ 
systems with — 1372 for several cooling rates. We 
find that Gg (r) plateaus at large r and the plateau value 
G'g(oo) increases with decreasing cooling rate R. For 
partially crystalhne systems, Gg(r) decays to I/VTV^ at 
large distances, where is the number of independent 
crystalline domains. For disordered systems, Gg(r) de- 
cays to where Nj, is the total number of near- 
est neighbor particles [Ij]. We find that the deviation 
Gg(rniax) — Gg(oo), where Gg(rinax) are the local maxima 
in Gg(r), decays exponentially ^ e~^^^ with correlation 
length^. (See Fig. [5]) The correlation length ^ grows lin- 
early with the linear size of the system iV^/'^ for cooling 
rates R < Rc- 

We also employed a crystal analysis algorithm to iden- 
tify the crystalhne clusters (FCC, HCP [IJI, or BCC) 
that form during the thermal quenching process. For 
slow cooling rates, the system forms only a few large 
crystalline clusters whose size scales with the system size. 
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FIG. 8: (color online) (left) Local bond orientational order correlation function Ge{r) for monodisperse LJ systems with 
— 1372 at several cooling rates -R = 1, 10~^, 10~^, and 10~^. (middle) The decay of the local maxima in G6(r) versus 
distance r for monodisperse LJ systems at cooling rate R = 10"'^ for several system sizes A'" = 500, 864, 1372, 2048, and 4000. 
(right) Correlation length ^ from the decay of the local bond orientational order correlation function versus the linear dimension 
of the system A'^^''^ for monodisperse LJ systems at cooling rate R = 10~^. The solid line has slope ^ 0.5. 



(See Fig. [5]). For fast cooling rates, the number of crys- 
talline clusters is small, and each cluster contains only 
a few particles. At intermediate rates, the number of 
crystalline clusters reaches a maximum at a characteris- 
tic cooling rate that scales with N as shown in Fig. [SI 



These results are consistent with the fact that the criti- 
cal cooling rate Rc (defined using the local bond orienta- 
tional order parameter Qg) becomes independent of the 
system size in the N ^ oo limit. 



[1] G. Kumar, P. Neibecker, Y. H. Liu, and J. Schroers, Na- 
ture Communications 4, 1536 (2013). 
[2] A. Inoue, Acta Mater. 48, 279 (2000). 
[3] D. TurnbuU, Contemp. Phys. 10, 473 (1969). 
[4] A. Takeuchi and A. Inoue, Mater. Trans. 46, 2817 (2005). 
[5] D. B. Miracle, W. S. Sanders, and O. N. Senkov, Philos. 

Mag. 83, 2409 (2003). 
[6] P. Jalah and M. Li, Intermetallics 12, 1167 (2004). 
[7] P. Jalah and M. Li, Phys. Rev. B 71, 014206 (2005). 
[8] M. P. Allen and D.J. Tildesley, Computer Simulation of 

Liquids (Oxford University Press, New York, 1987). 
[9] A. B. Hopkins, Y. Jiao, F. H. Stillinger, and S. Torquato, 
Phy. Rev. Lett. 107, 125501 (2011). 
[10] S. Nose, J. Ghem. Phys. 81, 511 (1984). 
[11] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). 
[12] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. 

Rev. B 28, 784 (1983). 
[13] C. F. Schreck and G. S. O'Hern, in Experimental and 
Computational Techniques in Soft Condensed Matter 
Physics, edited by J. Olafsen (Cambridge University 
Press, Cambridge, 2010), pp. 25-61. 
[14] Y. T. Wang, S. Teitel, and C. Dellago, J. Ghem. Phys. 

122, 214722 (2005). 
[15] A. Stukowski, Modelhng Simul. Mater. Sci. Eng. 20, 
045021 (2012). 

[16] HCP-like particles and clusters reported in this work are 
identified by considering the first nearest neighbors of 
each particle. If instead, the super-lattice across stacking 
layers was included, most of the HCP clusters would be 
classified as '9R' structures |29j] (with repeating motifs 



composed of an FGC layer followed by two HGP layers). 
[17] M. D. Rintoul and S. Torquato, Phys. Rev. Lett. 77, 4198 
(1996). 

[18] J. D. Honeycutt and H. C. Andersen, Ghem. Phys. Lett. 

108, 535 (1984). 
[19] J. D. Honeycutt and H. C. Andersen, J. Phys. Ghem. 90, 

1585 (1986). 

[20] A. R. Miedema, R. Boom, and F. R. Deboer, J. Less- 

Gommon Met. 41, 283 (1975). 
[21] R. Boom, F. R. de Boer, and A. R. Miedema, J. Less- 

Gommon Met. 45, 237 (1976). 
[22] R. Boom, F. R. de Boer, and A. R. Miedema, J. Less- 

Gommon Met. 46, 271 (1976). 
[23] J. H. O. Varley, Philos. Mag. 45, 887 (1954). 
[24] Y. Q. Cheng, E. Ma, and H. W. Sheng, Phys. Rev. Lett. 

102 (2009). 

[25] X. J. Liu, X. D. Hui, G. L. Chen, and T. Liu, Phys. Lett. 

A 373, 2488 (2009). 
[26] O. N. Senkov, Y. Q. Cheng, D. B. Miracle, E. R. Barney, 

A. C. Hannon, and C. F. Woodward, J. App. Phys. Ill, 

123515 (2012). 

[27] D. Frenkel and B. Smit, Understanding Molecular Simu- 
lation (Academic Press, 2002). 

[28] D. Brown and J. H. R. Clarke, Mol. Phys. 51, 1243 
(1984). 

[29] F. Ernst, M. W. Finnis, D. Hofmann, T. Muschik, 
U. Schonberger, and U. Wolf, Phys. Rev. Lett. 69, 620 
(1992). 



8 




lO'* 10"^ 10"^ 10' 10° lO'* 10^ 10"^ 10"' 10° 

R R 



FIG. 9: (color online) (left) The number of crystalline clusters Uc (FCC, HCP, and BCC) normalized by the system size A*' 
for monodisperse LJ systems as a function of cooling rate R for several system sizes, (right) The number of (FCC, HCP, and 
BCC) crystal-like particles Nc normalized by the number of crystalline clusters ric {i.e. average crystalline cluster size) as a 
function of cooling rate for several system sizes. 



